home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Turnbull China Bikeride
/
Turnbull China Bikeride - Disc 1.iso
/
ARGONET
/
PD
/
MATHS
/
RLAB
/
RLAB125.ZIP
/
!RLaB
/
toolbox
/
section
< prev
next >
Wrap
Text File
|
1994-09-23
|
4KB
|
129 lines
//-------------------------------------------------------------------//
// Synopsis: Calculate the section properties (area, center of
// geometry, moment of inertia) of a region.
// Syntax: section ( P , coord )
// Description:
// The properties of polygonal sections may be calculated using this
// rfile.
//
// The (x,y) coordinates of the vertices of the region are stored
// in P. P is a two column matrix, its first column stores xs and
// its second column stores ys. The coordinates of the vertces must
// be stored sequentially for a complete clockwise path around the
// region. Holes inside the region can be subtracted by tracing
// around them in counter clockwise direction.
// If coord is "polar", then the coordinates in P are in polar
// coordinates (r,theta), where theta is in degree.
//
// Reference: (many more)
// Wojciechowski, F., "Properties of Plane Cross Sections,"
// Machine Design, Jan. 1976, p.105.
//
// Example 1: Section properties of a rectangle
//
// A = [0,0;0,4;2,4;2,0;0,0];
// S = section(A);
// show_prop(A,S);
// pause();
//
// Example 2: Section properties of a triangle
//
// P = [0,0;0,1;1,0;0,0];
// S = section(P);
// show_prop(P,S);
// pause();
//
// Example 3: Section properties of an unit circle
//
// theta = (360:0:-5)';
// r = ones(theta.nr,1);
// S = section([r,theta],"polar");
// show_prop([r,theta],S,"polar");
// pause();
//
// Example 4: Section properties of a C section
//
// P = [1,75;ones(42,1),(80:285:5)';1.2,285;ones(42,1)*1.2,(280:75:-5)'];
// S = section(P,"polar");
// show_prop(P,S,"polar");
// pause();
//
// Example 5: Section properties of a square with a center hole
// s2 = sqrt(2);
// P = [10,0;10*s2,-45;10*s2,225;10*s2,135;10*s2,45;10,0;...
// 5*ones(73,1),(0:360:5)';10,0];
// S = section(P,"polar");
// show_prop(P,S,"polar");
// pause();
//
// Todo:
// 1. make P a linked-list of several regions (may be holes), return
// the total properties.
// 2. add thickness and mass density to each region in 1., so the mass,
// the center of mass, and the mass moment of inertia can be
// calculated.
//
// Tzong-Shuoh Yang 8/10/94 (tsyang@ce.berkeley.edu)
//-------------------------------------------------------------------//
section = function( P, coord )
{
local(P,B,N,N1,np,prop,q);
global (eps, pi)
if (P.nr <= 2)
{
error("section.r: A polygon must have at least 3 points.");
}
// if the polygon is not closed, let the last point = 1st point
if (P[P.nr;1] != P[1;1] || P[P.nr;2] != P[1;2] )
{
P = [P;P[1;1],P[1;2]];
}
// convert polar coordinates to cartesian coordinates
if (exist(coord))
{ if (coord == "polar")
{
P[;2] = P[;2]*(pi/180);
P = [P[;1].*cos(P[;2]),P[;1].*sin(P[;2])];
}
}
np = P.nr - 1;
N = 1:np;
N1 = 2:P.nr;
prop = <<>>;
B = P[N1;1].*P[N;2] - P[N;1].*P[N1;2];
prop.area = sum(B)/2;
prop.cgx = B'*(P[N;1]+P[N1;1])/6/prop.area;
prop.cgy = B'*(P[N;2]+P[N1;2])/6/prop.area;
prop.Ixx = B'*(P[N;2].^2 + P[N;2].*P[N1;2] + P[N1;2].^2)/12;
prop.Iyy = B'*(P[N;1].^2 + P[N;1].*P[N1;1] + P[N1;1].^2)/12;
prop.Ixy = B'*(P[N;1].*(2*P[N;2]+P[N1;2]) + P[N1;1].*(2*P[N1;2]+P[N;2]))/24;
// centroid values
prop.Ixx0 = prop.Ixx - prop.cgy^2*prop.area;
prop.Iyy0 = prop.Iyy - prop.cgx^2*prop.area;
prop.Ixy0 = prop.Ixy - prop.cgx*prop.cgy*prop.area;
// pricipal values
if (abs(prop.Ixy0) < eps) {
q = 0;
else if (prop.Ixx0 == prop.Iyy0) {
q = pi/4;
else
q = 0.5*atan(2*prop.Ixy0/(prop.Iyy0 - prop.Ixx0));
}}
prop.pIxx = prop.Ixx0*cos(q)^2 + prop.Iyy0*sin(q)^2 - prop.Ixy0*sin(2*q);
prop.pIyy = prop.Ixx0*sin(q)^2 + prop.Iyy0*cos(q)^2 + prop.Ixy0*sin(2*q);
// polar moi
prop.J0 = prop.pIxx + prop.pIyy;
// convert angle to degree
prop.angle = q*180/pi;
return prop;
};